rm(list=ls())

pacman::p_load(lubridate, tidyverse, haven, 
               gtrendsR, rdrobust, rddtools,
               broom, ggpubr, gridExtra,readxl)

setwd('put_your_wd_here')

allshootings <- read_csv('datasets/list_of_shootings.csv',col_select = 1:4)
allshootings <- allshootings %>%
  filter(include == 1)
top10 <- allshootings$shooting[allshootings$top10 == 1]

# this function was used to pull gtrends data from the 
# google trends API, estimate RDITs for each shooting
# and save the results

# gtrends_table <- function(term, date, shooting, time=2){
#   
#   cut = lubridate::ymd(date) 
#   start_date = cut - months(time)
#   end_date = cut + months(time)
#   time_span = paste(start_date, end_date, sep=' ')
#   res <- gtrends(term,geo='US',
#                   low_search_volume = T,
#                   time = time_span)
#   dat <- res$interest_over_time
#   dat$date <- ymd(dat$date)
#   dat$running_day <- as.numeric(dat$date - cut)
#   dat$hits <- as.numeric(dat$hits)
#   dat <- na.omit(dat)
#   
#   sd_est <- sd(as.numeric(dat$hits),na.rm=T)
#   house_rdd <- rdd_data(y=as.numeric(dat$hits), x=dat$running_day, cutpoint=0)
#   bw_ik <- rdd_bw_ik(house_rdd)
#   reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
#   rddtoolsest <- data.frame(Est = 'Rdd Tools Coef',
#                            Coef=reg_nonpara$coefficients/sd_est,
#                            lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2])/sd_est,
#                            upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2])/sd_est,
#                            shooting=shooting,
#                            term=term)
# 
#   out <- bind_rows(rddtoolsest)
#   return(out)
# }

# this code uses the function from above and pulls data

# terms <- c('Everytown for Gun Safety','Brady Campaign',
#            'Gun Control','Gun Rights',
#            'NRA','Gun Owners of America')
# 
# results <- vector('list', nrow(allshootings)*length(terms))
# count = 0
# for(j in 1:nrow(allshootings)){
#   for(i in 1:length(terms)){
#     out <- tryCatch({gtrends_table(term=terms[i],
#                                   date=allshootings$date[j],
#                                   shooting=allshootings$shooting[j])},
#                     error = function(e){NA})
#     out
#     count = count + 1
#     if(is.na(out)){
#       results[[count]] <- NA
#     } else {
#     results[[count]] <- out
#     filename <- paste0(allshootings$shooting[j],'_',terms[i],'.csv')
#     write_csv(dat, file = glue::glue('/Users/tyler/Dropbox/Mass Shootings and Tweets, Trends, and Donations/data/trends/raw_data2/{filename}'))
#     }
#     print(paste0("shooting ",j,", term ",i, ", count ",count))
#   }
# }
# 
# dat <- bind_rows(results[lengths(results) > 1])
# write_csv(dat, 'trends_dat_full_new.csv')
# 
# # placebos 
# terms <- c('Abortion',
#            'Recycling',
#            'Sexism')
# 
# results <- vector('list', nrow(allshootings)*length(terms))
# count = 0
# for(j in 1:nrow(allshootings)){
#   for(i in 1:length(terms)){
#     out <- tryCatch({gtrends_table(term=terms[i],
#                                    date=allshootings$date[j],
#                                    shooting=allshootings$shooting[j])},
#                     error = function(e){NA})
#     out
#     count = count + 1
#     if(is.na(out)){
#       results[[count]] <- NA
#     } else {
#       results[[count]] <- out
#       filename <- paste0(allshootings$shooting[j],'_',terms[i],'.csv')
#       write_csv(dat, file = glue::glue('{filename}'))
#     }
#     print(paste0("shooting ",j,", term ",i, ", count ",count))
#   }
# }
# dat <- bind_rows(results[lengths(results) > 1])
# write_csv(dat, 'trends_dat_full_new_placebos.csv')

dat <- read_csv('datasets/trends_dat_full_allshootings.csv')
dat$shooting[dat$shooting == 'Aurora Hentry Pratt'] <- 'Aurora Henry Pratt'

# split into groups
dat_ideas <- dat %>% filter(term %in% c('Gun Control','Gun Rights'))
dat_gc <- dat %>% filter(term %in% c('Everytown for Gun Safety','Brady Campaign'))
dat_gr <- dat %>% filter(term %in% c('NRA','Gun Owners of America'))
placebos <- dat %>% filter(term  == 'Recycling')

# lower 20 meta for figures
lower20 <- read_csv('datasets/lower_20_google.csv')
lower20$shooting[lower20$shooting=='Meta Estimate (All)'] <- 'Meta Estimate (Bottom 20)'

# meta analysis function
calcMeta <- function(data, name, shooting, term){
  m.gen <- meta::metagen(TE = Coef,
                         seTE = se,
                         studlab = shooting,
                         data = data,
                         sm = "SMD",
                         fixed = FALSE,
                         random = TRUE)
  results <- tibble(shooting = name, 
                    Coef = m.gen$TE.random,
                    lower = m.gen$lower.random,
                    upper = m.gen$upper.random,
                    term = term)
  return(results)
}

# append meta analyses
top20 <- allshootings$shooting[allshootings$include==1]
meta1 <- calcMeta(dat_ideas %>% 
                   filter(term == 'Gun Rights' & shooting %in% top20) %>%
                   mutate(se = (upper - Coef)/1.96), 
                 'Meta Estimate (Top 20)','Meta Analysis', 'Meta Gun Rights')
meta2 <- calcMeta(dat_ideas %>% 
                    filter(term == 'Gun Control' & shooting %in% top20) %>%
                    mutate(se = (upper - Coef)/1.96), 
                  'Meta Estimate (Top 20)','Meta Analysis', 'Meta Gun Control')
meta3 <- calcMeta(dat_ideas %>% 
                    filter(term == 'Gun Rights' & shooting %in% top10) %>%
                    mutate(se = (upper - Coef)/1.96), 
                  'Meta Estimate (Top 10)','Meta Analysis', 'Meta Gun Rights')
meta4 <- calcMeta(dat_ideas %>% 
                    filter(term == 'Gun Control' & shooting %in% top10) %>%
                    mutate(se = (upper - Coef)/1.96), 
                  'Meta Estimate (Top 10)','Meta Analysis', 'Meta Gun Control')
dat_ideas <- bind_rows(dat_ideas, meta1, meta2, meta3, meta4, lower20[1:2,])

# save meta analytic effects
# bind_rows(meta1, meta2, meta3, meta4, lower20[1:2,]) %>%
#  mutate(z = 'Gun Control and Gun Rights Searches Google Trends') %>%
#  write_csv(file='datasets/metaanalysis/meta_google_all.csv')

# Gun Control Orgs
meta3a <- calcMeta(dat_gc %>% 
                    filter(shooting %in% top20) %>%
                    mutate(se = (upper - Coef)/1.96), 
                  'Meta Estimate (Top 20)','Meta Analysis', 'Meta')
meta3b <- calcMeta(dat_gc %>% 
                    mutate(se = (upper - Coef)/1.96) %>%
                     filter(shooting %in% top10), 
                  'Meta Estimate (Top 10)','Meta Analysis', 'Meta')
dat_gc <- bind_rows(dat_gc, meta3a, meta3b,lower20[3,])

# bind_rows(meta3a, meta3b,lower20[3,]) %>%
#   mutate(z = 'Gun Control Org Searches Google Trends') %>%
#   write_csv(file='datasets/metaanalysis/meta_google_gun_control_orgs.csv')

# Ideas
meta4a <- calcMeta(dat_gr %>% 
                     filter(shooting %in% top20) %>%
                    mutate(se = (upper - Coef)/1.96), 
                  'Meta Estimate (Top 20)','Meta Analysis', 'Meta')
meta4b <- calcMeta(dat_gr %>% 
                    mutate(se = (upper - Coef)/1.96) %>%
                     filter(shooting %in% top10), 
                  'Meta Estimate (Top 10)','Meta Analysis', 'Meta')
dat_gr <- bind_rows(dat_gr, meta4a, meta4b, lower20[4,])

# bind_rows(meta4a, meta4b,lower20[4,]) %>%
#   mutate(z = 'Gun Rights Org Searches Google Trends') %>%
#   write_csv(file='datasets/metaanalysis/meta_google_gun_rights_orgs.csv')

####
# Meta analyses by race of victims
#### 

race_est <- read_csv('datasets/shooting_race.csv')
race_est <- race_est %>%
  filter(shooting %in% dat$shooting)
terciles <- quantile(race_est$mean.white, na.rm=T, c(0,0.33,0.66,1))
race_est$tercile <- ifelse(race_est$mean.white < terciles[2],1,0)
race_est$tercile <- ifelse(race_est$mean.white > terciles[2] & race_est$mean.white < terciles[3],2,race_est$tercile)
race_est$tercile <- ifelse(race_est$mean.white > terciles[3],3,race_est$tercile)
datnew <- left_join(dat, race_est, by='shooting')

meta_top3 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 3 & term %in% c('Gun Owners of America','NRA')),  
                      'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis')
meta_top2 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 2 & term %in% c('Gun Owners of America','NRA')),  
                      'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis')
meta_top1 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 1  & term %in% c('Gun Owners of America','NRA')),  
                      'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis')
meta_tercile_gunrightsorgs <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
  mutate(category = 'Google', z='Gun Rights Orgs')

meta_top3 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 3 & term %in% c('Everytown for Gun Safety','Brady Campaign')),  
                      'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis')
meta_top2 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 2 & term %in% c('Everytown for Gun Safety','Brady Campaign')),  
                      'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis')
meta_top1 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 1  & term %in% c('Everytown for Gun Safety','Brady Campaign')),  
                      'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis')
meta_tercile_guncontrolorgs <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
  mutate(category = 'Google', z='Gun Control Orgs')

meta_top3 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 3 & term == 'Gun Control'),  
                      'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis')
meta_top2 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 2 & term == 'Gun Control'), 
                      'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis')
meta_top1 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 1  & term == 'Gun Control'),
                      'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis')
meta_tercile_guncontrol <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
  mutate(category = 'Google', z='Gun Control')

meta_top3 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 3 & term == 'Gun Rights'),  
                      'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis')
meta_top2 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 2 & term == 'Gun Rights'), 
                      'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis')
meta_top1 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 1  & term == 'Gun Rights'),
                      'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis')
meta_tercile_gunrights <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
  mutate(category = 'Google', z='Gun Rights')
meta_tercile_out <- bind_rows(
  meta_tercile_gunrightsorgs,
  meta_tercile_guncontrolorgs,
  meta_tercile_guncontrol,
  meta_tercile_gunrights
)
  
write_csv(meta_tercile_out, 
          'datasets/metaanalysis/meta_race_google.csv' )

# making plots
dat_ideas_c <- dat_ideas %>% 
   mutate(shooting = case_when(
     shooting == "Aurora Movie Theater" ~ "Aurora Movie Theater (2012)",
     shooting == "Sandy Hook" ~ "Sandy Hook (2012)",
     shooting == "Washington Navy Yard" ~ "Washington Navy Yard (2013)",
     shooting == "Isla Vista" ~ "Isla Vista (2014)",
     shooting == "Charleston Church" ~ "Charleston Church (2015)",
     shooting == "Umpqua" ~ "Umpqua (2015)",
     shooting == "San Bernardino" ~ "San Bernardino (2015)",
     shooting == "Orlando" ~ "Orlando (2016)",
     shooting == "Las Vegas" ~ "Las Vegas (2017)",
     shooting == "Sutherland Springs Church" ~ "Sutherland Springs Church (2017)",
     shooting == "Stoneman Douglas" ~ "Stoneman Douglas (2018)",
     shooting == "Santa Fe High" ~ "Santa Fe High (2018)",
     shooting == "Pittsburgh Synagogue" ~ "Pittsburgh Synagogue (2018)",
     shooting == "Thousand Oaks Borderline Bar" ~ "Thousand Oaks Borderline Bar (2018)",
     shooting == "Aurora Henry Pratt" ~ "Aurora Henry Pratt (2019)",
     shooting == "Virginia Beach Municipal" ~ "Virginia Beach Municipal (2019)",
     shooting == "Gilroy + El Paso + Dayton" ~ "Gilroy + El Paso + Dayton (2019)",
     shooting == "Atlanta + Boulder" ~ "Atlanta + Boulder (2021)",
     shooting == "Indianapolis Fed Ex" ~ "Indianapolis Fed Ex (2021)",
     shooting == "San Jose VTA" ~ "San Jose VTA (2021)",
     TRUE ~ shooting
   ), shooting = factor(shooting, levels=rev(c(
     "Aurora Movie Theater (2012)",
       "Sandy Hook (2012)",
       "Washington Navy Yard (2013)",
       "Isla Vista (2014)",
       "Charleston Church (2015)",
       "Umpqua (2015)",
       "San Bernardino (2015)",
       "Orlando (2016)",
       "Las Vegas (2017)",
       "Sutherland Springs Church (2017)",
       "Stoneman Douglas (2018)",
       "Santa Fe High (2018)",
       "Pittsburgh Synagogue (2018)",
       "Thousand Oaks Borderline Bar (2018)",
       "Aurora Henry Pratt (2019)",
       "Virginia Beach Municipal (2019)",
       "Gilroy + El Paso + Dayton (2019)",
       "Atlanta + Boulder (2021)",
       "Indianapolis Fed Ex (2021)",
       "San Jose VTA (2021)",
      "Meta Estimate (Top 10)",
     "Meta Estimate (Top 20)",
     "Meta Estimate (Bottom 20)")))) 

dat_ideas_c$black <- ifelse(
  (dat_ideas_c$lower < 0 & dat_ideas_c$upper < 0) |
    (dat_ideas_c$lower > 0 & dat_ideas_c$upper > 0),1,0)

plot_ideas <- dat_ideas_c %>%
  mutate(term = factor(term, levels=c(
    'Gun Rights','Gun Control','Meta Gun Rights','Meta Gun Control'
  ))) %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, shape=term, color=factor(black))) +
  geom_errorbar(width=0, position=position_dodge(width=0.5)) +
  geom_point(position=position_dodge(width=0.5), size=2, fill='white') +
  coord_flip() +
  geom_hline(yintercept=0, linetype=2, color='red') +
  labs(title='(B) Policy Positions', x='',y='') +
  scale_color_manual(values=c('grey50','Black')) +
  scale_shape_manual(values=c(24, 21, 17, 16)) +
  theme_bw() +
  theme(legend.position='none', 
        plot.title = element_text(hjust = 0.5))  
plot_ideas

# plot 2 placebos
meta5a <- calcMeta(placebos %>% 
                     filter(shooting %in% allshootings$shooting[allshootings$include==1]) %>%
                     mutate(se = (upper - Coef)/1.96),
                   'Meta Estimate (Top 20)','Meta Analysis', 'Meta')
meta5b <- calcMeta(placebos %>% 
                     filter(shooting %in% allshootings$shooting[allshootings$include==1] &
                              shooting %in% top10) %>%
                     mutate(se = (upper - Coef)/1.96),
                   'Meta Estimate (Top 10)','Meta Analysis', 'Meta')
placebos <- bind_rows(placebos, meta5a, meta5b, lower20[5,])

missing_shoots <- allshootings$shooting[!allshootings$shooting %in% placebos$shooting]
placebos <- bind_rows(placebos, tibble(shooting = missing_shoots)) %>%
  filter(shooting %in% c(allshootings$shooting,'Meta Estimate (Top 10)','Meta Estimate (Top 20)','Meta Estimate (Bottom 20)'))
placebos <- placebos %>% 
  mutate(shooting = factor(shooting, levels=rev(c(
    "Aurora Movie Theater",
    "Sandy Hook",
    "Washington Navy Yard",
    "Isla Vista",
    "Charleston Church",
    "Umpqua",
    "San Bernardino",
    "Orlando",
    "Las Vegas",
    "Sutherland Springs Church",
    "Stoneman Douglas",
    "Santa Fe High",
    "Pittsburgh Synagogue",
    "Thousand Oaks Borderline Bar",
    "Aurora Henry Pratt",
    "Virginia Beach Municipal",
    "Gilroy + El Paso + Dayton",
    "Atlanta + Boulder",
    "Indianapolis Fed Ex",
    "San Jose VTA",
    "Meta Estimate (Top 10)",
    "Meta Estimate (Top 20)",
    "Meta Estimate (Bottom 20)")))) 

placebos$black <- ifelse(
  (placebos$lower < 0 & placebos$upper < 0) |
    (placebos$lower > 0 & placebos$upper > 0),1,0)

plot_pl <- placebos %>%
  filter(!is.na(shooting)) %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, shape=term,color=factor(black))) +
  coord_flip() +
  geom_hline(yintercept=0, linetype=2, color='red') +
  labs(title='(C) Placebo', x='',y='') +
  scale_color_manual(values=c('grey50','Black')) +
  theme_bw() +
  theme(legend.position='none', 
        plot.title = element_text(hjust = 0.5),
        axis.text.y = element_blank()) +
  geom_errorbar(width=0, position=position_dodge(width=0.25)) +
  scale_shape_manual(values=c(19, 21)) +
  geom_point(position=position_dodge(width=0.25), size=2,  fill='white') 
plot_pl

dat_gc <- dat_gc %>% 
  mutate(shooting = case_when(
  shooting == "Aurora Movie Theater" ~ "Aurora Movie Theater (2012)",
  shooting == "Sandy Hook" ~ "Sandy Hook (2012)",
  shooting == "Washington Navy Yard" ~ "Washington Navy Yard (2013)",
  shooting == "Isla Vista" ~ "Isla Vista (2014)",
  shooting == "Charleston Church" ~ "Charleston Church (2015)",
  shooting == "Umpqua" ~ "Umpqua (2015)",
  shooting == "San Bernardino" ~ "San Bernardino (2015)",
  shooting == "Orlando" ~ "Orlando (2016)",
  shooting == "Las Vegas" ~ "Las Vegas (2017)",
  shooting == "Sutherland Springs Church" ~ "Sutherland Springs Church (2017)",
  shooting == "Stoneman Douglas" ~ "Stoneman Douglas (2018)",
  shooting == "Santa Fe High" ~ "Santa Fe High (2018)",
  shooting == "Pittsburgh Synagogue" ~ "Pittsburgh Synagogue (2018)",
  shooting == "Thousand Oaks Borderline Bar" ~ "Thousand Oaks Borderline Bar (2018)",
  shooting == "Aurora Henry Pratt" ~ "Aurora Henry Pratt (2019)",
  shooting == "Virginia Beach Municipal" ~ "Virginia Beach Municipal (2019)",
  shooting == "Gilroy + El Paso + Dayton" ~ "Gilroy + El Paso + Dayton (2019)",
  shooting == "Atlanta + Boulder" ~ "Atlanta + Boulder (2021)",
  shooting == "Indianapolis Fed Ex" ~ "Indianapolis Fed Ex (2021)",
  shooting == "San Jose VTA" ~ "San Jose VTA (2021)",
  TRUE ~ as.character(shooting)
),shooting = factor(shooting, levels=rev(c(
  "Aurora Movie Theater (2012)",
  "Sandy Hook (2012)",
  "Washington Navy Yard (2013)",
  "Isla Vista (2014)",
  "Charleston Church (2015)",
  "Umpqua (2015)",
  "San Bernardino (2015)",
  "Orlando (2016)",
  "Las Vegas (2017)",
  "Sutherland Springs Church (2017)",
  "Stoneman Douglas (2018)",
  "Santa Fe High (2018)",
  "Pittsburgh Synagogue (2018)",
  "Thousand Oaks Borderline Bar (2018)",
  "Aurora Henry Pratt (2019)",
  "Virginia Beach Municipal (2019)",
  "Gilroy + El Paso + Dayton (2019)",
  "Atlanta + Boulder (2021)",
  "Indianapolis Fed Ex (2021)",
  "San Jose VTA (2021)",
  "Meta Estimate (Top 10)",
  "Meta Estimate (Top 20)",
  "Meta Estimate (Bottom 20)"))))

dat_gc$black <- ifelse(
  (dat_gc$lower < 0 & dat_gc$upper < 0) |
    (dat_gc$lower > 0 & dat_gc$upper > 0),1,0)

plot_gc <- dat_gc %>%
   ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, shape = term, color=factor(black))) +
   geom_errorbar(width=0, position=position_dodge(width=0.5)) +
   geom_point(position=position_dodge(width=0.5), size=2, fill='white') +
   coord_flip() +
   geom_hline(yintercept=0, linetype=2, color='red') +
   labs(title='(D) Gun Control Orgs', x='',y='') +
   scale_color_manual(values=c('grey50','Black')) +
   scale_shape_manual(values=c(24,21,16)) +
   theme_bw() +
   theme(legend.position='none', 
         plot.title = element_text(hjust = 0.5)) 
plot_gc

# gun rights
 dat_gr <- dat_gr %>% 
   mutate(shooting = factor(shooting, levels=rev(c(
     "Aurora Movie Theater",
     "Sandy Hook",
     "Washington Navy Yard",
     "Isla Vista",
     "Charleston Church",
     "Umpqua",
     "San Bernardino",
     "Orlando",
     "Las Vegas",
     "Sutherland Springs Church",
     "Stoneman Douglas",
     "Santa Fe High",
     "Pittsburgh Synagogue",
     "Thousand Oaks Borderline Bar",
     "Aurora Henry Pratt",
     "Virginia Beach Municipal",
     "Gilroy + El Paso + Dayton",
     "Atlanta + Boulder",
     "Indianapolis Fed Ex",
     "San Jose VTA",
     "Meta Estimate (Top 10)",
     "Meta Estimate (Top 20)",
     "Meta Estimate (Bottom 20)")))) 
 
 dat_gr$black <- ifelse(
   (dat_gr$lower < 0 & dat_gr$upper < 0) |
     (dat_gr$lower > 0 & dat_gr$upper > 0),1,0)
 
 plot_gr <- dat_gr %>%
   mutate(
     term = factor(term, levels=c('NRA','Gun Owners of America','Meta')),
          ) %>%
   ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, shape=term, color=factor(black))) +
   coord_flip() +
   geom_hline(yintercept=0, linetype=2, color='red') +
   labs(title='(E) Gun Rights Orgs', x='',y='') +
   scale_color_manual(values=c('grey50','Black'),guide='none') +
   theme_bw() +
   geom_errorbar(width=0, position=position_dodge(width=0.5)) +
   geom_point(position=position_dodge(width=0.5), size=2,fill='white') +
   scale_shape_manual(values=c(24,21,16),name='') + 
  theme(legend.position='none', 
        plot.title = element_text(hjust = 0.5),
        axis.text.y = element_blank())
plot_gr

plots <- ggarrange(plot_ideas, plot_pl,plot_gc,plot_gr,
                    heights = c(4,4,4,4),
                    widths = c(9,5,9,5),
                    ncol = 2, nrow = 2)
plots
plots <- annotate_figure(plots, bottom = text_grob("Effect on Google Searches\n(Standard Deviations)",family='Helvetica'))
ggsave(plots, width=8, height=8, 
       file='figures/gsearch_all2_new.png', bg='white')
 
#############
# bottom 20 results plots -----
#############

rm(list=ls())

allshootings <- read_csv('datasets/list_of_shootings.csv',col_select = 1:4)
allshootings <- allshootings %>%
  filter(include == 0)
allshootings$date <- mdy(allshootings$date)
allshootings$shooting[allshootings$shooting == 'Lafayette Move Theater'] <- "Lafayette Movie Theater"

# low attn plots
dat <- read_csv('datasets/trends_dat_full_lowattshootings.csv')
dat$shooting[dat$shooting == 'Lafayette Move Theater'] <- 'Lafayette Movie Theater'
dat$term <- sapply(str_split(dat$term, '_'), function(x) x[2])
dat$term <- gsub('.csv','',dat$term)

# subset
append_missing <- function(x){
  terms = unique(x$term)
  n = length(terms)
  list.out = vector('list', length=n)
  missing_shoots <- unique(allshootings$shooting[!allshootings$shooting %in% x$shooting])
  for(i in 1:n){
    list.out[[i]] <- tibble(Est = as.character(NA), Coef = NA, lower = NA, upper = NA, shooting = missing_shoots, term = terms[i])
  }
  results <- bind_rows(list.out)
  return(results)
}

dat_ideas <- dat %>% filter(term %in% c('Gun Control','Gun Rights')) 
dat_ideas <- dat_ideas %>%
  bind_rows(append_missing(dat_ideas))
dat_gc <- dat %>% filter(term %in% c('Everytown for Gun Safety','Brady Campaign'))
dat_gc <- dat_gc %>%
  bind_rows(append_missing(dat_gc))
dat_gr <- dat %>% filter(term %in% c('NRA','Gun Owners of America'))
dat_gr <- dat_gr %>%
  bind_rows(append_missing(dat_gr))
placebos <- dat %>% filter(term  == 'Recycling')
placebos <- placebos %>%
  bind_rows(append_missing(placebos))

# meta analysis
calcMeta <- function(data, name, shooting, term){
  m.gen <- meta::metagen(TE = Coef,
                         seTE = se,
                         studlab = shooting,
                         data = data,
                         sm = "SMD",
                         fixed = FALSE,
                         random = TRUE)
  results <- tibble(shooting = name, 
                    Coef = m.gen$TE.random,
                    lower = m.gen$lower.random,
                    upper = m.gen$upper.random,
                    term = term)
  return(results)
}

# Ideas
meta1 <- calcMeta(dat_ideas %>% 
                    filter(term == 'Gun Rights' ) %>%
                    mutate(se = (upper - Coef)/1.96), 
                  'Meta Estimate (All)','Meta Analysis', 'Meta Gun Rights')
meta2 <- calcMeta(dat_ideas %>% 
                    filter(term == 'Gun Control') %>%
                    mutate(se = (upper - Coef)/1.96), 
                  'Meta Estimate (All)','Meta Analysis', 'Meta Gun Control')
dat_ideas <- bind_rows(dat_ideas, meta1, meta2)

# Gun Control Orgs
meta3a <- calcMeta(dat_gc %>% 
                     mutate(se = (upper - Coef)/1.96), 
                   'Meta Estimate (All)','Meta Analysis', 'Meta')
dat_gc <- bind_rows(dat_gc, meta3a)

# Ideas
meta4a <- calcMeta(dat_gr %>% 
                     mutate(se = (upper - Coef)/1.96), 
                   'Meta Estimate (All)','Meta Analysis', 'Meta')
dat_gr <- bind_rows(dat_gr, meta4a)

# plot 1
dat_ideas_c <- dat_ideas %>%
  mutate(shooting = factor(shooting, levels=rev(c(
    allshootings$shooting, 'Meta Estimate (All)'
  ))))

dat_ideas_c$black <- ifelse(
  (dat_ideas_c$lower < 0 & dat_ideas_c$upper < 0) |
    (dat_ideas_c$lower > 0 & dat_ideas_c$upper > 0),1,0)

plot_ideas <- dat_ideas_c %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, shape=term, color=factor(black))) +
  geom_errorbar(width=0, position=position_dodge(width=0.5)) +
  geom_point(position=position_dodge(width=0.5), size=2, fill='white') +
  coord_flip() +
  geom_hline(yintercept=0, linetype=2, color='red') +
  labs(title='(B) Policy Positions', x='',y='') +
  scale_color_manual(values=c('grey50','Black')) +
  scale_shape_manual(values=c(24, 21, 17, 16)) +
  theme_bw() +
  theme(legend.position='none', 
        plot.title = element_text(hjust = 0.5))  
plot_ideas

# plot 2
dat_gc <- dat_gc %>% 
  mutate(shooting = factor(shooting, levels=rev(c(
    allshootings$shooting,"Meta Estimate (All)")))) 

dat_gc$black <- ifelse(
  (dat_gc$lower < 0 & dat_gc$upper < 0) |
    (dat_gc$lower > 0 & dat_gc$upper > 0),1,0)

plot_gc <- dat_gc %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, shape = term, color=factor(black))) +
  geom_errorbar(width=0, position=position_dodge(width=0.5)) +
  geom_point(position=position_dodge(width=0.5), size=2, fill='white') +
  coord_flip() +
  geom_hline(yintercept=0, linetype=2, color='red') +
  labs(title='(C) Gun Control Orgs', x='',y='') +
  scale_color_manual(values=c('grey50','Black')) +
  scale_shape_manual(values=c(24,21,16)) +
  theme_bw() +
  theme(legend.position='none', 
        plot.title = element_text(hjust = 0.5),
        axis.text.y = element_blank())
plot_gc

# gun rights
dat_gr <- dat_gr %>% 
  mutate(shooting = factor(shooting, levels=rev(c(
    allshootings$shooting,"Meta Estimate (All)")))) 

dat_gr$black <- ifelse(
  (dat_gr$lower < 0 & dat_gr$upper < 0) |
    (dat_gr$lower > 0 & dat_gr$upper > 0),1,0)

plot_gr <- dat_gr %>%
  mutate(
    term = factor(term, levels=c('NRA','Gun Owners of America','Meta')),
  ) %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, shape=term, color=factor(black))) +
  coord_flip() +
  geom_hline(yintercept=0, linetype=2, color='red') +
  labs(title='(D) Gun Rights Orgs', x='',y='') +
  scale_color_manual(values=c('grey50','Black'),guide='none') +
  theme_bw() +
  theme(legend.position='none', 
        plot.title = element_text(hjust = 0.5),
        axis.text.y = element_blank()) +
  geom_errorbar(width=0, position=position_dodge(width=0.5)) +
  geom_point(position=position_dodge(width=0.5), size=2,fill='white') +
  scale_shape_manual(values=c(24,21,16),name='')
plot_gr

#  placebos
meta5a <- calcMeta(placebos %>% 
                     filter(shooting %in% allshootings$shooting[allshootings$include==0]) %>%
                     mutate(se = (upper - Coef)/1.96),
                   'Meta Estimate (All)','Meta Analysis', 'Meta')
placebos <- bind_rows(placebos, meta5a)

# write out meta analysis lower 20 to file
# write_csv(bind_rows(meta1, meta2, meta3a, meta4a,meta5a) %>%
#             mutate(term2=c('Gun Rights','Gun Control',
#                            'Gun Control Orgs','Gun Rights Orgs',
#                            'Placebo')),
#           file = 'datasets/metaanalysis/lower_20_google.csv')

placebos <- placebos %>% 
  mutate(shooting = factor(shooting, levels=rev(c(
    allshootings$shooting,"Meta Estimate (All)")))) 

placebos$black <- ifelse(
  (placebos$lower < 0 & placebos$upper < 0) |
    (placebos$lower > 0 & placebos$upper > 0),1,0)

plot_pl <- placebos %>%
  filter(!is.na(shooting)) %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, shape=term,color=factor(black))) +
  coord_flip() +
  geom_hline(yintercept=0, linetype=2, color='red') +
  labs(title='(E) Placebo', x='',y='') +
  scale_color_manual(values=c('grey50','Black')) +
  theme_bw() +
  theme(legend.position='none', 
        plot.title = element_text(hjust = 0.5),
        axis.text.y = element_blank()) +
  geom_errorbar(width=0, position=position_dodge(width=0.25)) +
  scale_shape_manual(values=c(19, 21)) +
  geom_point(position=position_dodge(width=0.25), size=2,  fill='white') 
plot_pl

plots <- ggarrange(plot_ideas, plot_gc,plot_gr,plot_pl,
                   heights = c(4,4,4,4),
                   widths = c(9,5,5,5),
                   ncol = 4, nrow = 1)
plots
plots <- annotate_figure(plots, bottom = text_grob("Effect on Google Searches\n(Standard Deviations)",family='Helvetica'))
ggsave(plots, width=12, height=5, file='gsearch_all2_low_attn.png', bg='white')
